#include <time.h>
#include <sys/time.h>
#include <math.h>
#include "../include/mersenneTwisterCPU.h"
#include "../include/rw_functions.h"
#include "../include/dci.h"

# define DIV 16

static mt_struct 	 MT[MT_RNG_COUNT];

static uint32_t state[32][MT_NN*16];

/*
 * Função:   rw_serial
 * Entradas: Número de Particulas, Dimensão da matriz, Semente para geração de 
             números aleatórios.
 * Objetivo: Executar o algoritmo de RW de forma Serial
 * Saida:    Uma matriz com todo o calculo executado
 */
int* rw_serial( const int NParticles, 
                const int Size,
			    const int seed )
{
    int *mat;
    mat = (int*) calloc (Size * Size, sizeof(int) );    

    //Contole da poscicao da particula no dominio
    int x, y;
    int random_number, ok;
    
    //variaveis de contole
    int particles, step;
    int parar = Size*Size;
    //double start, end;
    const int RNG_COUNT = NParticles / DIV;
    int iRng;
    
    
    if ( (NParticles % DIV) != 0 ){
        printf("Atenção o número de particulas deve ser maior que 16 e potencia de 2\n");
        printf("Número de Particulas Atual: %d", NParticles ); 
        printf("Ex: 16, 32, 64, 128, 256, 512, 1024...\n");
        exit(EXIT_FAILURE);
    }
    
    
    //start = omp_get_wtime();
    
    for ( iRng = 0; iRng < RNG_COUNT; iRng++ )
    {
        MT[iRng].state = state;
        sgenrand_mt(seed, &MT[iRng]);
        
        for( particles = 0; particles < DIV; particles++ ) 
        {
            /* colocando a particula no ponto inicial */
            x = y = 0;

            ++mat[0];
            step = 0;
                    
            /*
             * iniciando a caminhada de uma particula
             * Parada por distancia percorrida ou por chegada em casa
             */
            
            while(step < parar)
            {
                random_number = genrand_mt(&MT[iRng]) % 4;

                ok = 0;

                switch ( random_number ){
                    case 0:
                        if ( ( x - 1 ) >= 0 ){
                            --x;
                            ok = 1;
                        }
                        break;
                    case 1:
                        if ( ( y + 1 ) < Size ){
                            ++y;
                            ok = 1;
                        }
                        break;
                    case 2:
                        if ( ( x + 1) < Size ){
                            ++x;
                            ok = 1;
                        }
                        break;
                    case 3:
                        if ( (y -1) >= 0 ){
                            --y;
                            ok = 1;
                        }
                        break;
                }
                if(ok){
                    step++;
                    ++mat[(x * Size) + y];
                    
                }

				if ((x==Size-1)&&(y==Size-1))
					break;
            }
        }
    }
    //end = omp_get_wtime();
    //printf("%lf\t%f\n", end-start, (mat[Size*Size-1] / (double)NParticles) );

    return mat;
}


/*
 * Função:   rw_pOpenMP
 * Entradas: Número de Particulas, Dimensão da matriz, Semente para geração de 
             números aleatórios, número de threads de Execução 
 * Objetivo: Executar o algoritmo de RW de forma paralela com OpenMP
 * Saida:    Uma matriz com todo o calculo executado
 */
int* rw_pOpenMP( const int NParticles, 
                 const int Size, 
                 const unsigned int seed, 
                 const int numThreads )
{
    int *mat;
    mat = (int*) calloc (Size * Size, sizeof(int) ); 
    
    //Contole da poscicao da particula no dominio
    int x, y;
    //Controle da poscicao na matriz
    int pos;
    //variaveis de contole
    int particles, step;
    int iRng;
    int random_number, ok;
    int parar = Size*Size;
    const int RNG_COUNT = NParticles / DIV;
    int my_thread_id;
    
    if ( (NParticles % DIV) != 0 ){
        printf("Atenção o número de particulas deve ser um valor potência de 2\n");
        printf("Número de Particulas Atual: %d", NParticles ); 
        printf("Ex: 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024...\n");
        exit(EXIT_FAILURE);
    }
    
    omp_set_num_threads(numThreads);
    
    double start, end;
                 
    start = omp_get_wtime(); /* Inicializa captação do tempo */
    
    #pragma omp parallel private(my_thread_id)
	{
		my_thread_id = omp_get_thread_num();
		        
		#pragma omp for private(iRng, particles, step, pos, x, y, random_number, ok)
		for ( iRng = 0; iRng < RNG_COUNT; iRng++ )
		{                
		    MT[iRng].state = &state[my_thread_id][0];
		    sgenrand_mt(seed, &MT[iRng]);
		    
		    for( particles = 0; particles < DIV; particles++ ) 
		    {   
		        x = y = 0;
		        
		        #pragma omp atomic
		        mat[0]++;
		        
		        step = 0;

		        while(step < parar)
		        {
		            random_number = genrand_mt(&MT[iRng]) % 4; //Gera números Aleatórios Apartir do código do Mersenne Twister
		            ok = 0;

		            switch ( random_number ){
		                case 0:
		                    if ( ( x - 1 ) >= 0 ){
		                        --x;
		                        ok = 1;
		                    }
		                    break;
		                case 1:
		                    if ( ( y + 1 ) < Size ){
		                        ++y;
		                        ok = 1;
		                    }
		                    break;
		                case 2:
		                    if ( ( x + 1) < Size ){
		                        ++x;
		                        ok = 1;
		                    }
		                    break;
		                case 3:
		                    if ( (y -1) >= 0 ){
		                        --y;
		                        ok = 1;
		                    }
		                    break;
		            }
		            if(ok){
		                step++;
				        pos = (x * Size) + y;
		                
		                #pragma omp atomic
		                mat[pos]++;
		                
		                if (pos == (parar-1) ) 
		                    break;
		            }
		        }
		    }  
		}
    }
    
    end = omp_get_wtime(); /* Finaliza captação do tempo */
    printf("%lf\t%f\n", end-start, (mat[Size*Size-1] / (double)NParticles) );
    return mat;
}

/*
 * Função: rw_pOpenMPCrit
 * Entradas: Número de Particulas, Dimensão da matriz, Semente para geração de 
             números aleatórios, número de threads de Execução 
 * Objetivo: Executar o algoritmo de RW de forma paralela com OpenMP
 * Saida:    Uma matriz com todo o calculo executado
 */
int* rw_pOpenMPCrit( const int NParticles, 
                 const int Size, 
                 const unsigned int seed, 
                 const int numThreads )
{
    int *mat;
    mat = (int*) calloc (Size * Size, sizeof(int) ); 
    
    //Contole da poscicao da particula no dominio
    int x, y;
    //Controle da poscicao na matriz
    int pos;
    //variaveis de contole
    int iParticulas, step;
    int iRng;
    int random_number, ok;
    int parar = Size*Size;
    const int RNG_COUNT = NParticles / DIV;
    
    
    if ( (NParticles % DIV) != 0 ){
        printf("Atenção o número de particulas deve ser um valor potência de 2\n");
        printf("Número de Particulas Atual: %d", NParticles ); 
        printf("Ex: 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024...\n");
        exit(EXIT_FAILURE);
    }
    
    omp_set_num_threads(numThreads);
    
    double start, end;
                 
    start = omp_get_wtime(); /* Inicializa captação do tempo */
            
    #pragma omp parallel for private(iRng, iParticulas, step, pos, x, y, random_number, ok)
    for ( iRng = 0; iRng < RNG_COUNT; iRng++ )
    {                
        MT[iRng].state = &state[omp_get_thread_num()][0];
        sgenrand_mt(seed, &MT[iRng]);
        
        for( iParticulas = 0; iParticulas < DIV; iParticulas++ ) 
        {   
            x = y = 0;
            
            #pragma omp critical
            ++mat[0];
            
            step = 0;

            while(step < parar)
            {
                random_number = genrand_mt(&MT[iRng]) % 4; //Gera números Aleatórios Apartir do código do Mersenne Twister
                ok = 0;

                switch ( random_number ){
                    case 0:
                        if ( ( x - 1 ) >= 0 ){
                            --x;
                            ok = 1;
                        }
                        break;
                    case 1:
                        if ( ( y + 1 ) < Size ){
                            ++y;
                            ok = 1;
                        }
                        break;
                    case 2:
                        if ( ( x + 1) < Size ){
                            ++x;
                            ok = 1;
                        }
                        break;
                    case 3:
                        if ( (y -1) >= 0 ){
                            --y;
                            ok = 1;
                        }
                        break;
                }
                if(ok){
                    step++;
		            pos = (x * Size) + y;
		            
                    #pragma omp critical
                    ++mat[pos];
                    
                    if ( pos == (parar-1) ) 
                        break;
                }
            }
        }  
    }
    
    end = omp_get_wtime(); /* Finaliza captação do tempo */
    printf("Tempo Paralelo: %lf\nProbabilidade: %f\n", end-start, (mat[Size*Size-1] / (double)NParticles) );
    return mat;
}


/*
 * Função: rw_parallel_reduction
 * Entradas: Número de Particulas, Dimensão da matriz, Semente para geração de 
             números aleatórios, número de threads de Execução 
             
 * Objetivo: Executar o algoritmo de RW de forma paralela com OpenMP utilizando mais de uma matriz
             a fim de evitar o uso de operações atomicas. Obs. O número de matrizes alocadas é
             diretamente proporcional ao número de Threads disponíveis.
             Após a execussão do RW, é feita uma redução de todas as matrizes para apenas uma
             matriz que é retornada
             
 * Saida:    Uma matriz com todo o calculo executado
 */
int* rw_parallel_reduction( const int NParticles, 
                            const int Size, 
                            const unsigned int seed, 
                            const int numThreads )
{
    int jCount, iCount;
    int **mat = NULL;
    mat = (int**) calloc ( numThreads, sizeof(int*) );
    
    for ( iCount = 0; iCount < numThreads; iCount++ )
        mat[iCount] = (int*) calloc (Size * Size, sizeof(int) );    

    int *matriz = NULL;
    matriz = (int*) calloc (Size * Size, sizeof(int) ); 
    
    omp_set_num_threads(numThreads);

    //Contole da poscicao da particula no dominio
    int x, y, iRng;
    int random_number, ok;
    
    //variaveis de contole
    int particles, step;
    int parar = Size*Size;
    const int RNG_COUNT = NParticles / DIV;
    
    if ( (NParticles % DIV) != 0 ){
        printf("Atenção o número de particulas deve ser um valor potência de 2\n");
        printf("Número de Particulas Atual: %d", NParticles ); 
        printf("Ex: 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024...\n");
        exit(EXIT_FAILURE);
    }
    
    int threadAtual;
         
    //double start, end;
                 
    //printf("Paralelo OpenMp Reduction, Número de Particulas: %d\n", NParticles );
    //start = omp_get_wtime();
    
    #pragma omp parallel for private(iRng, particles, step, x, y, random_number, ok)
    for ( iRng = 0; iRng < RNG_COUNT; iRng++ )
    {
        //int prova = 0;
        threadAtual = omp_get_thread_num();
        MT[iRng].state = &state[threadAtual][0];
        sgenrand_mt(seed, &MT[iRng]);
            
        for( particles = 0; particles < DIV; particles++ ) 
        {
            //prova++;

            /* colocando a particula no ponto inicial */
            x = y = 0;

            ++mat[threadAtual][0];
            step = 0;
                    
            /*
             * iniciando a caminhada de uma particula
             * Parada por distancia percorrida ou por chegada em casa
             */
            
            while(step < parar)
            {
                random_number = 1 + genrand_mt(&MT[iRng]) % 4;

                ok = 0;

                switch ( random_number ){
                    case 1:
                        if ( ( x - 1 ) >= 0 ){
                            --x;
                            ok = 1;
                        }
                        break;
                    case 2:
                        if ( ( y + 1 ) < Size ){
                            ++y;
                            ok = 1;
                        }
                        break;
                    case 3:
                        if ( ( x + 1) < Size ){
                            ++x;
                            ok = 1;
                        }
                        break;
                    case 4:
                        if ( (y -1) >= 0 ){
                            --y;
                            ok = 1;
                        }
                        break;
                }
                if(ok){
                    step++;
                    ++mat[threadAtual][(x * Size) + y];
                    
                    if ((x == (Size - 1))&&(y == (Size - 1))) 
                        break; 
                }
            }
        }
        
        //printf("Fim Thread: %d\tProva: %d\n", threadAtual, prova );
    
    } /* Fim da região paralela */
    
    
    #pragma omp parallel for private (jCount, iCount)
    for ( iCount = 0; iCount < Size * Size; iCount++ )
        for ( jCount = 0; jCount < numThreads; jCount++ )
            matriz[iCount] += mat[jCount][iCount];
    
    
    //end = omp_get_wtime();
    //printf("Chegaram: %d\nNumParticulas %d\n", matriz[Size*Size-1], NParticles );
    //printf("%lf\t%lf\n", end-start, (matriz[Size*Size-1] / (double)NParticles) );

    return matriz;
}


/**************************************************************************************************/
/* Métodos iterativos que tem como objetivo encontrar a propabilidade exata para cada esperimento */


/*

 * Função:   rw_serial_iterativo
 * Entradas: taxa de incremento, Dimensão da matriz, Semente para geração de 
             números aleatórios 
 * Objetivo: Executar o algoritmo de RW de forma serial e iterativo afim de obter o número perfeito de
             particulas que devem ser disparadas 

 * Saida:    Uma matriz com todo o calculo executado
 */

int* rw_serial_iterativo ( int incremento, 
                const int Size,
			    int seed )
{
    int *mat;
    mat = (int*) calloc (Size * Size, sizeof(int) );    

    //Contole da poscicao da particula no dominio
    int x, y;
    int random_number, ok;
    
    //variaveis de contole
    int particles, step;
    int parar = Size*Size;
    double start, end, oldProb = 0.0, newProb = 0.0, error = 0.0;
    const int taxaIncremento = incremento;
    const int RNG_COUNT = taxaIncremento / DIV;
    int iRng, numInteration = 0, flag = 0;
    
    if ( (taxaIncremento % DIV) != 0 ){
        printf("Atenção a taxa de incremento deve ser um valor potência de 2\nTaxa de Incremento atual: %d\nEx: 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024...\n", incremento );
        exit(EXIT_FAILURE);
    }
    
    start = omp_get_wtime();
    
    do
    { 
        for ( iRng = 0; iRng < RNG_COUNT; iRng++ )
        {
            MT[iRng].state = state;
            sgenrand_mt(seed, &MT[iRng]);
            
            for( particles = 0; particles < DIV; particles++ ) 
            {
                /* colocando a particula no ponto inicial */
                x = y = 0;

                ++mat[0];
                step = 0;
                        
                /*

                 * iniciando a caminhada de uma particula
                 * Parada por distancia percorrida ou por chegada em casa
                 */
                
                while(step < parar)
                {
                    random_number = genrand_mt(&MT[iRng]) % 4;

                    ok = 0;

                    switch ( random_number ){
                        case 0:
                            if ( ( x - 1 ) >= 0 ){
                                --x;
                                ok = 1;
                            }
                            break;
                        case 1:
                            if ( ( y + 1 ) < Size ){
                                ++y;
                                ok = 1;
                            }
                            break;
                        case 2:
                            if ( ( x + 1) < Size ){
                                ++x;
                                ok = 1;
                            }
                            break;
                        case 3:
                            if ( (y -1) >= 0 ){
                                --y;
                                ok = 1;
                            }
                            break;
                    }
                    if(ok){
                        step++;
                        ++mat[(x * Size) + y];
                        
                        if ((x == (Size - 1))&&(y == (Size - 1))) 
                            break;
                    }
                }
            }
        }
        
        printf("Chegaram: %d\nNumParticulas %d\n", mat[Size*Size-1], incremento );
        newProb = (mat[Size*Size-1] / (double)incremento);
        
        error = newProb - oldProb;
        
        error = ( error > 0.000 ) ? error : (-1) * error;
        printf("Taxa de Erro: %lf\n", error );
        
        if ( error < (1e-3) )
            flag = 1;
        
        numInteration++;
        oldProb = newProb;
        
        seed = time(NULL);
        incremento = incremento + taxaIncremento;
        
        printf("Probabilidade: %lf\nNumero de Iterações: %d\n", newProb, numInteration ); 
        
    } while ( flag == 0 );
    
    end = omp_get_wtime();
    
    printf("Tempo Serial: %lf\n", end-start );
    return mat;
}


/*

 * Função:   rw_pOpenMP_iterativo
 * Entradas: Taxa de incremento, Dimensão da matriz, Semente para geração de 
             números aleatórios, número de Threads 
 * Objetivo: Executar o algoritmo de RW de forma paralela com OpenMp de forma iterativa afim de obter
  o número perfeito de particulas que devem ser disparadas 

 * Saida:    Uma matriz com todo o calculo executado
 */
int* rw_pOpenMP_iterativo ( int incremento, 
                const int Size,
			    unsigned int seed,
			    const int numThreads )
{
    int *mat;
    mat = (int*) calloc (Size * Size, sizeof(int) );    

    //Contole da poscicao da particula no dominio
    int x, y;
    int random_number, ok;
    
    //variaveis de contole
    int iParticulas, step, pos;
    int parar = Size*Size;
    double start, end, oldProb = 0.0, newProb = 0.0, error = 0.0;
    const int taxaIncremento = incremento;
    const int RNG_COUNT = taxaIncremento / DIV;
    int iRng, numInteration = 0, flag = 0;
    
    
    if ( (taxaIncremento % DIV) != 0 ){
        printf("Atenção a taxa de incremento deve ser um valor potência de 2\n");
        printf("Ex: 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024...\n");
        exit(EXIT_FAILURE);
    }
    
    omp_set_num_threads(numThreads);
                 
    start = omp_get_wtime(); /* Inicializa captação do tempo */
    
    do {
            
        #pragma omp parallel for private(iRng, iParticulas, step, pos, x, y, random_number, ok)
        for ( iRng = 0; iRng < RNG_COUNT; iRng++ )
        {                
            MT[iRng].state = &state[omp_get_thread_num()][0];
            sgenrand_mt(seed, &MT[iRng]);
            
            for( iParticulas = 0; iParticulas < DIV; iParticulas++ ) 
            {   
                x = y = 0;
                
                #pragma omp atomic
                mat[0]++;
                
                step = 0;

                while(step < parar)
                {
                    random_number = genrand_mt(&MT[iRng]) % 4; //Gera números Aleatórios Apartir do código do Mersenne Twister
                    ok = 0;

                    switch ( random_number ){
                        case 0:
                            if ( ( x - 1 ) >= 0 ){
                                --x;
                                ok = 1;
                            }
                            break;
                        case 1:
                            if ( ( y + 1 ) < Size ){
                                ++y;
                                ok = 1;
                            }
                            break;
                        case 2:
                            if ( ( x + 1) < Size ){
                                ++x;
                                ok = 1;
                            }
                            break;
                        case 3:
                            if ( (y -1) >= 0 ){
                                --y;
                                ok = 1;
                            }
                            break;
                    }
                    if(ok){
                        step++;
		                pos = (x * Size) + y;
                        #pragma omp atomic
                        mat[pos]++;
                        
                        if ((x == (Size - 1))&&(y == (Size - 1))) 
                            break;
                    }
                }
            }  
        }
    
        printf("Chegaram: %d\nNumParticulas %d\n", mat[Size*Size-1], incremento );
        newProb = (mat[Size*Size-1] / (double)incremento);
        
        error = newProb - oldProb;
        
        error = ( error > 0.000 ) ? error : (-1) * error;
        printf("Taxa de Erro: %lf\n", error );
        
        if ( error < (1e-3) )
            flag = 1;
        
        numInteration++;
        oldProb = newProb;
        
        seed = time(NULL);
        incremento = incremento + taxaIncremento;
        
        printf("Probabilidade: %lf\nNumero de Iterações: %d\n\n", newProb, numInteration ); 
            
    } while ( flag == 0 );
    
    end = omp_get_wtime();
    printf ("Tempo Paralelo: %lf\n", end - start );
    
    return mat;
}




void initMTRefCPU(const char *fname )
/* Initialize MT structures from file specified in 'fname'. If 'init_pgi_structures'
 * is TRUE, also initializes MT structures for PGI Accelerator directives. */
{
    FILE* fd = 0;
    if ((fd = fopen(fname, "rb")) == 0)
    {
        if(fd)
            fclose (fd);
    }

    for (int i = 0; i < MT_RNG_COUNT; i++) {
        //Inline structure size for compatibility,
        //since pointer types are 8-byte on 64-bit systems (unused *state variable)
        if( !fread(MT + i, 16 /* sizeof(mt_struct) */ * sizeof(int), 1, fd) ) {
//	        oclCheckError(0, 1);
        }
    }

    fclose(fd);
}
